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ABSTRACT 

Coalescing binary black holes experience a "kick" due to anisotropic emission 
of gravitational waves with an amplitude as great as ~ 200 km s _1 . We examine 
the orbital evolution of black holes that have been kicked from the centers of 
triaxial galaxies. Time scales for orbital decay are generally longer in triaxial 
galaxies than in equivalent spherical galaxies, since a kicked black hole does not 
return directly through the dense center where the dynamical friction force is 
highest. We evaluate this effect by constructing self-consistent triaxial models 
and integrating the trajectories of massive particles after they are ejected from 
the center; the dynamical friction force is computed directly from the velocity 
dispersion tensor of the self-consistent model. We find return times that are 
several times longer than in a spherical galaxy with the same radial density 
profile, particularly in galaxy models with dense centers, implying a substantially 
greater probability of finding an off-center black hole. 

Subject headings: black hole physics - galaxies: nuclei - galaxies: bulges - galaxies: 
kinematics and dynamics - galaxies: elliptical and lenticular, CD 



Introduction 



The coalescence of binary black holes (BHs) results in a gravitational recoil, or "kick" , 
due to the net linear momentum carried away by gravitational waves. Bekenstein (1973) 
estimated a kick velocity V ~ 300 km s _1 in highly nonspherical collapses using a quasi- 
Newtonian formalism. Nakamura et al. (1987) computed V/c — 0.045ry 2 for head-on col- 
lisions from infinity using BH perturbation theory; here rj = fi/(Mi + M 2 ), the "reduced 
mass ratio", with M\,Mi the BH masses and /i = \\\ \ / L >/( M\ + M2) the reduced mass. 
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A number of analytic estimates have been made of V for circular-orbit inspirals (Fitchett 
1983; Fitchett & Detweiler 1984; Wiseman 1992; Favata, Hughes & Holz 2004). The kick 
amplitude is usually divided into two components: the net recoil up to the innermost stable 
circular orbit (ISCO), and the contribution from the final plunge, from the ISCO to the 
horizon, which takes place in the strong-field regime and which dominates the total kick. 
Blanchet, Qusailah & Will (2005) computed V to second post-Newtonian order and found 
it to be well approximated by the simple formula 

^ = 0.04377V 1 - 4 *7 + (1) 

The i] 2 (l — At]) 1 / 2 dependence is the same found by Fitchett (1983), who computed gravita- 
tional recoil for a pair of BHs interacting via Newtonian forces and included only the lowest 
gravitational wave multipoles needed for momentum ejection. With the additional, ad hoc 
factor in equation (jTJ, Blanchet, Qusailah & Will (2005) were able to reproduce the results 
of their 2PN calculations to better than 1% at all mass ratios. The maximum estimated 
kick velocity was 250 ± 50 km s _1 . Damour & Gopakumar (2006) estimated a much lower 
value, ~ 74 km s _1 , for the maximum kick using an effective one-body approach. In the 
last year, remarkable progress has been made in techniques for the numerical solution of 
the full field equations (Pretorius 2005, 2006; Campanelli et al. 2005; Campanelli, Kelly & 
Lousto 2006; Baker et al. 2006a,b) allowing several groups to compute recoil velocities for 
coalescing BH binaries without approximations. Baker et al. (2006c) find V = 105 ± 10 km 
s _1 for M2/M1 = 2/3. Herrmann, Shoemaker & Laguna (2006) derive V = 33 km s _1 for 
M 2 /M l = 0.85 and V = 9 km s" 1 for M 2 /M l = 0.96. Most recently, Gonzalez et al. (2006) 
carried out a large set of inspiral simulations for non-spinning, circular-orbit binaries and 
determined the kick velocity as a function of mass ratio. Their results are well described by 
the expression 

V o / 

— = 0.040//V 1 - 477 (1 - 0.93r/) , (2) 

implying a maximum kick velocity of 175.7 km s _1 at M2/M1 = 0.36, somewhat smaller 
than implied by equation ([I]). The dependence of the kick velocity on orbital eccentricity 
was investigated by Sopuerta, Yunes and Laguna (2006) using the close-limit approximation, 
which models the late stages of coalescence as a single perturbed BH, coupled with a post- 
Newtonian estimate of the recoil during the early evolution. Their results are consistent 
with a (1 + e) dependence of kick velocity on eccentricity for the low (e < 0.1) eccentricities 
which they investigated, i.e. V ~ 195 km s _1 for e = 0.1. 

All of these results were obtained for non-spinning holes. In the presence of spin, kick 
velocities might be larger and would be nonzero even for Mi = M 2 (Redmount & Rees 
1989; Favata et al. 2004). Calculations of the coalescence of spinning BHs are currently 
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underway (Campanelli, Lousto and Zlochower 2006a, b). Recent results (Hermann et al. 
2007; Campanelli et al. 2007; Koppitz et al. 2007) actually show how recoil velocity after 
coalescence of spinning BHs may go up to ~ 450 km s -1 , reopening the possibility that a 
merged binary can be ejected even from the nucleus of a massive host galaxy. 

In this paper, we consider some of the consequences of the kicks for supermassive BHs 
in galactic nuclei. A kick velocity of 200 km s _1 is sufficient to remove a coalesced BH from 
a dwarf elliptical galaxy, even if the latter is embedded in a dark-matter halo (Merritt et al. 
2004). Escape velocities at the centers of luminous elliptical galaxies are generally greater 
than ~ 400 km s^ 1 however, and kicks of the magnitude so far calculated by numerical 
relativists would never be expected to remove BHs from such galaxies. But the kicks could 
still displace the BHs temporarily from their central locations, implying a finite probability 
of finding an off-center BH in a giant galaxy. The kicks would also generate long-lived 
changes in the central structure of galaxies as the displaced BHs transfer their orbital energy 
to the stars via dynamical friction (Merritt et al. 2004; Boylan-Kolchin et al. 2004; Madau 
& Quataert 2004). The displacements and their side-effects would have been greater at 
earlier times, when the gravitational potential wells associated with galaxies were shallower 
(Volonteri et al. 2003; Merritt et al. 2004; Madau & Quataert 2004; Haiman 2004; Yoo & 
Miralda-Escude 2004; Libeskind et al. 2006). 

In a non-spherical galaxy, an ejected BH does not pass precisely through the dense center 
as it falls back, reducing the mean value of the dynamical friction force compared with a 
spherical galaxy. The result is a more extended period of displacement compared with 
estimates based on spherical galaxy models. Here, we evaluate the effect of nonspherical 
galaxy geometries on BH infall times using fully self-consistent triaxial models. The models 
are constructed via orbital superposition as in Merritt & Fridman (1996), and the quantities 
that define the velocity ellipsoid at every point on the solution grid are computed and stored 
as in Merritt (1980). Given this information, it is possible to compute accurate estimates 
of the dynamical friction force that would act on a massive object, using the expressions 
in Pesce, Capuzzo-Dolcetta & Vietri (1992) and Capuzzo-Dolcetta & Vicari (2006). The 
frictional forces are added to the conservative forces from the triaxial mass distribution 
when integrating the BH trajectories. 

Among the new effects we find here is an increase in the effective "escape velocity" from 
the center of the galaxy, since some of the BH's initial kinetic energy is lost to dynamical 
friction. The kick velocity needed to escape is also dependent on the direction of the kick 
since the frictional force is direction-dependent. In general, we find that return times are 
longer in the triaxial geometry by factors of ~ a few compared with return times in spherical 
galaxies having the same, mean radial density profile, and this translates into a substantially 



-4- 



higher probability of finding a displaced BH. 



2. Method 

A BH ejected from the center of a galaxy moves in response to the conservative force 
corresponding to the gravitational potential $(r) of the stars and to the dynamical friction 
force per unit mass %(r). Its motion can be approximated by the solution of the differential 
equation 

r = -V$ + f df (3) 
subject to the proper initial conditions. We assume throughout that neither $ nor f^/ depend 
explicitly on time. 

As usual, we transform this second-order differential equation into a system of first-order 
equations: 

r = v, (4a) 
v = -V$ + %. (4b) 



As a consequence of the dynamical friction term, the orbital energy, E, is no longer a 
conserved quantity. Instead we can write 

E = E df = f df • v (5) 

This expression allows us to verify conservation of the total energy £ = K + $ + E d f, where 
K is the kinetic energy and Edf the work done by the dynamical friction force along the 
trajectory. 

We solved the system of differential equations (TjJ numerically, using the 7/8 order 
Runge-Kutta algorithm of Fehlberg (1968). □ 



2.1. Mass model 

We adopted as galaxy models the triaxial generalizations of the spherical Dehnen (1993) 
models. The mass density is 

(3-7)M 1 
Airaoc m 7 (l + m) 4 7 



1 The FORTRAN version of this algorithm was kindly made available to us by S. Udry. 
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where 

2 x 2 i y 2 ( z 2 



m^ = - + ^ + -, 0<c<6<a (7) 

and M is the total mass of the galaxy. The gravitational potential generated by this density 
law is (Chandrasekhar 1969) 



$(r) = -ixGabc 

Jo 



(oo)] — ^(m)]dr 



o ^/(r + a^^ + fo^^ + c 2 ) 



with 



and 



^(m) = / p(m' 2 )dm' 2 , (9) 
•0(oo) = lim ip(m), (10) 

m— *oo 

m 2 (r) = + -4- + (11) 

v 7 a 2 + r 6 2 + r c 2 + r 

Substituting s = yl + r in equation ([8]) leads (for 7 7^ 2) to an integral better suited to 
numerical evaluation (Merritt & Fridman 1996, hereafter MF96): 

= -^2 x 

x rl l-(3-7)( I jk) 2 -' + (2-7)( I fk) 3 ^ i „ fl2 X 

Jo v /[l+(6 2 -l)s 2 ][l+(c 2 -l)s 2 ] 1 ' 

with 

^ = S2[X2 + 1 + (J-1) 3 2 + l + (c ?- 1)s2 ]- ( 13 ) 

For 7 = 2 the potential is 

7 V /[1 + (62_1) S 2][ 1 + ( C 2_ 1 ) S 2] 

C = l09tdt (15) 

Jo v/[l + (6 2 -l)t 2 ][l + (c 2 -l)t 2 ]' 



The components of the force may be written 



(3 7) ~ ( ml {I + mi) 4 -V( a i + A* 2 )K 2 + ^s 2 ) (16) 
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where 



mf{s) 



x 



+ 



a 2 + Cis 2 a 2 + C 2 s 



+ 



C3S 2 



Here we have used the notation x± = x, x 2 = y, X3 
constants in equations ([TBI) and (IT7j) are 



Zj d\ — OS, Q2 



1 : 

Ai = b 2 - 1 A 2 
d = C 2 = b 2 



- 1 



(17) 

b, CI3 = c. The 



C 3 = c 2 - 1 



2 : 
Ai 



c 2 -6 2 
1-6 2 



A 2 = 1 - b 2 



Co = C^ = c 2 -b 2 



(18) 



3 : 
At 
Ci 



1-c 2 
1-c 2 



A 2 = 6 
C 2 = 6 2 



.2 ^2 



C, = 0. 



The integrals were computed by means of the double-precision FORTRAN routine DGAUSS 
of the CERNLIB library, which implements the classic Gaussian quadrature formula. 

Hereafter we adopt units in which G = M = a = 1. The unit of time is T u = a 3//2 / y/GM, 



or 



a 



3 / 2 / M \ -V2 



The unit of velocity is T4 = y/GM /a, or 



K = 666.8 A /—^ -kms -1 . 20 

a/kpc 



2.2. Self-consistent solutions 
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Fig. 1. — The radial behavior of the principal components of the velocity dispersion tensor 
in the three self-consistent models. The left column refers to model 1, the central column to 
model 2 and the right column to model 3 (see Tabled]). 
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The magnitude of the dynamical friction force depends on the details of the stellar 
velocity distribution. In order to accurately follow the orbital evolution of an ejected BH, a 
self-consistent galaxy model is therefore required. 

We constructed self-consistent triaxial models via a modification of the technique de- 
scribed in MF96. The densities generated by a catalogue of M orbits were recorded in a 
spatial grid of N = 1008 cells; the grid was defined as in MF96. We then sought the linear 
combination of orbits, with non-negative weights, which best reproduced the known mass in 
each cell imposing, as additional constraint, a zero streaming velocity. The quantity to be 
minimized was 



In this expression, Cj is the (unknown) number of stars on orbit i (1 < i < M); Bn is the 
mass which the ith orbit places in the /th cell; and D\ is the (known) mass which the model 
places in the Zth cell. The constraints Ci > were also imposed, i.e., the number of stars on 
each orbit was required to be positive. We used the NAG quadratic-programming routine 
E04NCF to carry out the minimization. 

We solved the self-consistency problem for the case 7 = 1, a "weak" density cusp, 
and 7 = 2, a "strong" cusp. In the 7=1 case, we built two models with different axis 
ratios, as reported in Table [TJ All models have a "triaxiality index" T = 0.5, where T = 
(a 2 — b 2 )/(a 2 — c 2 ); i.e. they are "maximally triaxial". The weak-cusp model is similar in 
structure to bright elliptical galaxies and bulges, those having absolute visual magnitudes 
brighter than My rs —21, while the strong-cusp model is similar to low- luminosity spheroids 
like that of the Milky Way and M32, which have higher- density nuclei. Table [D also gives, 
as useful reference time, the crossing time Ti/ 2 , defined as the period of the circular orbit at 
the radius containing one-half of the total galaxy mass, in the "equivalent" spherical model; 
the latter is defined as the spherical model with scale length given by (abc) 1 ^ 3 and the same 



In order to construct the velocity dispersion tensor, which is needed in the computation 
of the dynamical friction force, we also stored the velocities of the orbits as they passed 
through the cells. The velocity dispersion tensor is (Merritt 1980) 

N 




(21) 



1=1 



8=1 



total mass. 



J2^B a <V 2 k > 



u 




8=1 



(22) 



i=1 (^i-tSu 



where 



< Vjk >u=< VjVk >u ~ < Vj >u< V k > U ] 



(23) 
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here j, k refer to the coordinate axes, % is the orbit number, and I is the cell number. The 
principal components of the velocity dispersion tensor, erf, k = 1,2, 3, in each cell could then 
be computed via standard techniques. 

The high-energy cutoff in the orbital sample together with the unavoidable limitation 
in abundance of the orbital catalog implies that the outermost region of the spatial grid is 
not well populated, so we excluded from our evaluation of a^ k l the 48 cells of the outermost 
spatial shell. 

Figure [1] shows the radial behavior of the principal components of the velocity dispersion 
tensor in each of the three self-consistent models. 
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Table 1. 



Model 7 b/a c/a Ti/ 2 



1 1 0.7906 0.5 21.0 

2 1 0.8631 0.7 25.9 

3 2 0.8631 0.7 6.91 



Table 2. 



ID 


7 


r c /r B H a 


Core b 


v e c 


Ve,eff/V e d 


la 


1 





I 


1.625 


1.018 - 1.044 


Id 


1 


0.2 


II 


1.620 


1.008 - 1.036 


lc 


1 


2 


II 


1.605 


1.001 - 1.023 


If 


1 


0.2 


III 


1.620 


1.008 - 1.030 


lg 


1 


2 


III 


1.605 


1.001 - 1.017 


2a 


1 





I 


1.530 


1.030 - 1.047 


2d 


1 


0.2 


II 


1.530 


1.012 - 1.028 


2c 


1 


2 


II 


1.515 


1.004 - 1.019 


2f 


1 


0.2 


IV 


1.530 


1.009 - 1.024 


2g 


1 


2 


IV 


1.515 


1.004 - 1.013 


3b 


2 


0.2 


III 


4.66 


1.73 - 2.15 


3c 


2 


2 


III 


4.04 


1.89 - 2.24 


3d 


2 


0.2 


II 


4.66 


2.20 - 3.20 


3c 


2 


2 


II 


4.04 


1.11 - 1.35 


3f 


2 


0.2 


III 


4.66 


1.08 - 1.21 


3g 


2 


2 


III 


4.04 


1.00 - 1.05 



a Core radius in units of r BH (Eq. I2"9~j) . 

b Scheme for treatment of central forces (see §2.2). 

c Central escape velocity neglecting dynamical fric- 
tion. 

d Central escape velocity including dynamical fric- 
tion. 
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2.3. The dynamical friction force 

The classical Chandrasekhar (1943) formula for the dynamical friction force, which as- 
sumes a homogeneous and isotropic medium, has been extended by Pesce, Capuzzo-Dolcetta 
& Vietri (1992) to the triaxial case, in partial analogy with the Binney's treatment of the 
axysimmetric case (Binney 1977). The result is 

f dU = -T 1 e 1 - T 2 e 2 - r 3 e 3 (24) 

where (k = 1,2,3) are the principal components of the velocity dispersion tensor and 
Tfc(r, v, Vk) = 7fc(r, v)t> with Vk the component of the BH's velocity along the e^ axis. The 
functions 7,(r, v) are given by 

7i(r, v) = 2v / 27rp(r)G 2 In A(m + m,)a^ 3 x 

s 3 v k/ 2 "t 

r°° e e fc+" d (25) 

Here, m* is the mass of a field star, InA is the usual Coulomb logarithm, m. is the BH 
mass, G is the gravitational constant, p(r) is the galaxy mass density, and is the ratio 
between <7k and the largest eigenvalue, defined to be 03. 
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(x- (f) ,r (f) ) 

Fig. 2. — A two-dimensional illustration of the definition of x and y; the solid curve is the 
orbit of the BH. The ellipse has equation m 2 = m 2 . 
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The motion of the test object (the black hole) can be strongly affected by the values 
of the conservative and dynamical friction forces very near the center. Both components of 
the force would be influenced by the motion of the test particle; for instance, displacing a 
BH causes the central density to drop, an effect that was ignored in the construction of the 
models. In the case of the dynamical friction force, the limited resolution of our triaxial 
models is a problem; the central cells have a size of 0.28 and 0.05 for 7 = 1 and 7 = 2, 
respectively. We adopted the following scheme to deal with these issues. 

The standard Chandasekhar (1943) formula in the so called local approximation evalu- 
ates the dynamical friction force using the values of the stellar density and velocity at the 
position of the test particle. In the case of a particle that is at the center of a galaxy with 
a steep central density profile, the local approximation clearly yields an overestimate of the 
actual dynamical friction force, because it does not weight properly the contributions of par- 
ticles in zones surrouding the center, where the density is much smaller. To overcome this 
problem, Capuzzo-Dolcetta & Vicari (2006) proposed, for the case of a centrally-located test 
particle of mass m., the use of a numerical evaluation of the dynamical friction integral in 
its complete form. The fact that the test particle is at the center is exploited by identifying 
the distance to the scatterer (r) with the impact parameter b, yielding 

f df = -4vr ■ x 

m m + m 



Jv 1 + G 2 (m.+m) 2 



-G(v.)— , (26) 
v. 



where /(r, v) is the distribution function of the background particles (of mass m) that provide 
the frictional force. The integral in fT26|) is numerically performed assuming b m i n = and, 
for b max , a value large enough to guarantee convergence. 

After fitting with a suitable analytic expression the numerical G(v m ), we obtain an 
estimate of the dynamical friction force in the inner (m < m = 0.05) region of the galaxy 
by means of the interpolating formula 

fdf = -p{m)G{v.)— + [1 - p(m)]f df j(x, y, z), (27) 

where p(m) is a regular weighting function, satisfying p(0) = 1 and p(m) = 0. The expression 
[271 connects the central dynamical friction force, equation ff26l) . with its values at m > m, 
obtained via equation (j24"|) . Here, x(t),y(t) and z(t) are the coordinates of the intersection 
point between the radial direction through the position of the test particle at time t and the 
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m? = m? ellipsoidal surface (Figure [2]). In the following, we assumed for simplicity a linear 
dependence of p(m) on m. 

So far, we have ignored the influence of the BH and its motion on the distribution of 
mass in the galaxy. In reality, sudden removal of the BH from its central location will cause 
the stellar density to drop. The ejected BH will carry with it stars initially contained within 
a region r < r e // such that the orbital velocity around the BH is equal to Vkick, or 



the BH gravitational influence radius. Removal of the BH also reduces the gravitational 
force that binds the stars at r > tbh, causing them to move outward. Both effects result in 
a lowering of the density at r < r e ff. 

A self-consistent treatment of this expansion would require iV-body techniques. Here, 
we account approximately for the expansion by considering the dependence of our results 
on different assumptions about the inner form of the stellar density profile. We introduce 
a "core" radius r c , parametrized in terms of Tbh as r c = arsn, ot = (0.2,2). In computing 
r BH, w e set in Eq. [29] as a the value of the velocity dispersion in the innermost model cell. 




(28) 



with 



tbh = Gm./a 2 



(29) 
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Fig. 3. — Evolution of three planar (x = 0) orbits in model 2a (7 = 1), for V/v e = 
(0.4,0.5,0.6); for display convenience, we chose orbits with a different initial 6. Arrows 
indicate the direction of the motion for the two more energetic orbits. The lower panel 
shows the distance of the BH from the center as function of time. 
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Fig. 4. — As in Figure El for two planar orbits in model 3b (7 = 2), with V/v e = (0.8, 0.9). 
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Fig. 5. — Effective escape velocity as function of launching angle 9, for model la (left panel) 
and 3c (right panel). Empty circles refer to = 90°, crosses to = 45°, and black triangles 
to = 0°. 
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Note that tbh is 'under' the model resolution in all the cases treated, being tbh/L = 
0.16 and Tbh/L = 0.025 in the cases of 7 = 1 and 7 = 2, respectively (L is the innermost 
cell size, different in the two cases). 

Based on these considerations, we adopted four different schemes for treating the central 
forces (conservative + frictional), which we denote in Table [2] via the labels I-IV. These 
schemes are defined as follows. 

I No modification to the density profile, i.e. r c = 0. Conservative forces were computed 
as in equation ffTBI) (which is possible just when the central forces are finite, i.e. 7 = 
1) and dynamical friction forces as in equation (T2"T|) . i.e. by matching the triaxial 
expression to the central expression. 

II Conservative forces within the core were set to zero. Frictional forces in this region 
were computed via equation fl24|) . assuming for the coefficients their values at the 
boundary r = r c , as averaged over angles, and with the correct velocity dependence, 
i.e. T fc =< T k (r c ,v,v k ) ><^. 

III Density was set to zero within the core, thus conservative forces zero therein, and fric- 
tional forces in the core were evaluated setting b min = r c and p(m) = 1 in equation (T2"T|) . 

IV The density was given a core of radius r c and the potential within the core was assumed 
harmonic (linear forces), matching the external triaxial potential at the core boundary. 
Frictional forces were computed as in I. 

In what follows, we assumed a BH mass of m. = 10~ 3 in units of the galaxy mass, 
close to the mean ratio observed in real galaxies (Merritt & Ferrarese 2001; Marconi & Hunt 
2003). The Coulomb logarithm was set to In A = 6.6 (Spinnato et al. 2003). Dynamical 
friction times scale linearly with (m. InA) -1 . 

3. Results 

We studied the orbital evolution of the massive particle (BH) after being ejected from 
the center of each of the galaxy models in Table [2] with different kick velocities V and angles 
(0,4>), where V x = ^sin#cos0, V y = V sin 9 sin (j), V z = V^cos^. The kick velocity V was 
assigned a value in the range 0.2v e < V < v ej£ ff, where v e is the escape velocity from the 
origin in the absence of frictional forces, and f e , e // is the actual velocity required for escape; 
v e,eff > Ve since dynamical friction acts to reduce the particle's kinetic energy. Table [2] gives 
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v e and v e , e ff/v e (the latter expressed as a range, since f e , e // depends on the launching angle) 
for each of the models. The direction of the kick was given one of 43 values by choosing 9 
and (f) from the discrete set (0°, 15°, 30°, 45°, 60°, 75°, 90°). We defined the decay time T df as 
the time when the BH orbital energy had dropped to 1% its initial value and its residual 
time variation was negligible (\E/Eq\ < 10~ 4 ). 

Figures [3] and |4] show the evolution of a representative set of orbits. Figure [3] shows 
that for moderate kicks (V < 0.4t> e ) the BH executes only one "bounce" before dynamical 
friction brings it to a halt at the center. Figure 0] illustrates some higher energy orbits 
(V/v e = 0.8, 0.9) in a strong-cusp model (7 = 2). In spite of the strong frictional forces near 
the center of this model, ejection with a sufficiently large V allows the BH to execute several 
oscillations before coming to rest. 

The dependence of the effective escape velocity f e , e // 011 launching angle for two galaxy 
models is illustrated in Figure [51 The increase in u e ,e// relative to v e is most striking for 
trajectories aligned with the x— (long) axis, and in models with high central densities (7 = 2); 
v e,eff/ v e can be as large as ~ 3 (Table [2]). However, Table [2] also shows that v e>e ff/v e for 
these high density models can depend substantially on how the central forces are treated. 
For this reason, we consider the largest values of v e ^ e ff/v e in Table [2] to be provisional, until 
they can be verified with detailed iV-body simulations. 
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Fig. 6. — Left panel: Decay time T^f of the BH as function of the kick velocity V and for 
several launching angles (empty circles). In both panels, the dotted, short-dashed and long- 
dashed lines refer to kicks along the x-axis, y-axis and z-axis, respectively. The solid lines are 
for the "equivalent" spherical model, with scale length (a6c) 1//3 . Right panel: the maximum 
displacement of the BH. In both panels, dotted lines are for the "equivalent" spherical model, 
with scale length (abc) 1 ^ 3 . All results in this figure refer to model la. 
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Fig. 7. — Like Figure El but for model 2a. 




Fig. 9- 



Like Figure El but for model 3c 
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Fig. 10. — Orbital decay times vs the kick velocity for model If (dashed line) and lg (solid 
line) for the BH launched along the x (quickest decay), y and z (slowest decay) axis. 
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Fig. 11. — Orbital decay time in units of the period T of the unperturbed a;- axial orbit of 
the same energy, for models 2a (left panel) and 3b (right panel). 



10 3 



10 s 



-1 1 


1 1 
A 




i i i 


A 


i 

A 


1 






A 


□ 


U 


J 


A 




□ 


□ 






















J 


□ 








X 
























X 




O 


o 






X 


O 


o 






































O 
























O 


o 














"1 , 


I i 


1 I 


i i i 




1 


i 



0.5 1 




1.5 



10 3 



10 2 



-1 


6 


a 


i i 




i 


1 








H 












• 


A 










• 




O 


2 












• 














• 


□ 














A 
















A 














□ 












o 




• 










• 
















o- 














• 


"1 






1 1 




1 





0.5 1 1.5 



Fig. 12. — Orbital decay times in model la as a function of launching angle (left panel) 
and 9 (right panel), for vo/v e = 0.8. In each panel, symbols of the same type represent fixed 
values of the other angle, according to the scheme: black circles = 0°; black triangles = 30°; 
squares = 45°; crosses = 75°; empty circles=90°. 
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Figures [Ml] show the dependence of the decay time, Tdf, and of the apocentric "ellip- 
soidal" distance, m max , on the initial velocity in four of the models, for the entire set of 
launching angles. For the sake of comparison, the values of T d j and m max for the "equiva- 
lent" spherical model, characterized by the density law ([2D with r instead of m and length 
scale equal to (abc) 1 ^ 3 , are plotted as dotted lines. Kicks along the symmetry axes result 
in the shortest return times: the corresponding orbits are exactly radial and pass precisely 
through the dense center on each return, as in a spherical model. Return times are longer 
than in a spherical model, by as much as an order of magnitude in the 7 = 1 models, and 
even larger factors in the 7 = 2 models. In addition, Fig. [10] compares the decay times in 
two models (If and lg) with the same treatment of the dynamical friction term and same 
density profile but a different core size, a factor ten larger in mod. lg. Note the shorter 
decay times in the case of the smaller, higher density, core, due also to the smaller elongation 
of the radial oscillations at a given V/v e caused by the higher potential well. 

Figures 6-10 also show that oscillations along the major axis are quenched more rapidly 
than along the other demonstrated previously by Capuzzo-Dolcetta & Vicari (2005). 

Even though orbits launched close to the z- (short-) axis reach the greatest values of m max at 
every v/v e , their decay times are not necessarily the longest. This is because the dynamical 
friction force transforms the near-radial orbit into a box orbit which has a thin waist, thinner 
than that of orbits having the same initial energy but launching angles 9 7^ 0; the dynamical 
friction force is therefore greater because of the closer approach of the BH to the high-density 
regions. This is not a trivial result, because the galactic models treated in this paper are 
cuspy, and it was not a priori obvious that the large-scale behavior of the density would exert 
the dominant influence on the motion. At first sight it is surprising that the oscillations in the 
equivalent spherical models are always shorter than in the triaxial case, since in absence of 
dynamical friction, the radial oscillation in the spherical case reaches a maximum extension 
which is intermediate between the displacements along the x and z axes in the trixial case. 
However, the return time is always shortest in the spherical geometry since the effect of 
dynamical friction in the triaxial case is progressively greater from x to z, making the x 
oscillation the shortest, but still larger than in the spherical case. 

Decay times are often of the same order, or less than, the period of an unperturbed radial 
orbit in the spherical geometry, particularly in the models with 7 = 2. This is illustrated in 
Figure [11] for two models. Thus, in many cases, a kicked BH executes only ~ one or fewer 
complete radial oscillations before coming to rest. 

Figure [T2] illustrates the dependence of T d f on the launching angle. Decay times peak 
for intermediate angles, as was already apparent from Figures [ME 
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4. Consequences for Black Hole Displacements 

While triaxiality has the effect of lengthening the mean return time of a kicked BH to 
the center, compared with the time in a spherical galaxy (Figs. [Ml]), infall times are still 
typically short, of order a galaxy crossing time, unless the kick velocity is close to the escape 
velocity Another way to state this result is to say that there is a narrow range of kick 
velocities such that the BH spends a long time away from the center without being fully 
ejected. 

Here we estimate the probability that a kicked BH will be found at an appreciable 
distance from the center in a randomly-chosen galaxy. Since the distribution of kick velocities 
is poorly known, we will present results as a function of V/v e . The other parameter that 
determines the likelihood of finding a displaced BH is the time r since the galaxy experienced 
its last merger (which we equate with the elapsed time since the kick). This time is also 
poorly know for any galaxy. We therefore assume that r follows a Poisson distribution, 

g ^/^mcrgc 

p(r)dr = dr, (30) 

Emerge 

with t mer ge the mean time between mergers. 

Let -Py(r; At)dr be the probability of finding a kicked BH a distance r to r + dr from 
the center of the galaxy at a time At after the kick. In the case of a triaxial galaxy, we define 
Py as an average over the two launch angles (9, 0). Clearly for At > T d fy, this distribution 
is a delta function at the origin, where T^y is defined to be the maximum return time for 
kicks of magnitude V. In a randomly- chosen galaxy, the distribution of displacements for 
kicks of magnitude V is 

POO 

N v (r)= / P (T)P v (r;r)dT. (31) 



We simplify this expression by assuming that T d f is short compared with t merge , allowing 
us to approximate equation fl3~TT) by 

r T <if,v r°° 
N v (r)^P v (r) p(r)dT + 5(r) p(r)dr (32) 



where -FV( r ) is defined as the distribution of displacements averaged over the time < At < 
T dfy . Thus 

N V {r) W (l - e - T df,v/tmcr e ^ p y ^ + e -T d/ ,W*me Ig e£( r ) (33) 

and the cumulative distribution describing displacements less than r is 

N v (< r) » (1 - e - Td ^ v/tmcisc ) P v (< r) + e - T ^ /tmcrgc . (34) 
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We computed Py{< r) on a grid of V/v e - values for several of the models, then used 
these numbers to compute the probability distributions Ny(< r) as a function of the two 
parameters 

(Ti/ 2 /t mcTge ,V/v e ) (35) 

where Ty 2 is a proper reference time defined, as above (Tabled]), as the period of a circular 
orbit at the half-mass radius (~ T1/2) in the equivalent spherical model. Figures [TBI and [Til 
show the results for Models la and 2a; these models both have 7=1 and no "core" (Table [2]), 
differing only in their triaxiality (Table [TJ). These figures demonstrate that the probability 
of finding a kicked BH at a distance ~ ryi from the center of a galaxy is low unless the kick 
velocity is high, > 0.8w e , and the mean time between mergers is not too long compared with 
the crossing time. Based on Figs. [Mil this conclusion would be even stronger for models 
like 3b or 3c which have higher central densities. 
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Fig. 13. — Cumulative distributions iVy(< r) for the model la, for four different values of 
V/v e . Each frame plots equation ( 1341) for four different values of Ti/ 2 /t merge , where Ti/ 2 
is the period of a circular orbit at the half-mass radius in the equivalent spherical model, 
and irnerge ^ the mean time between galaxy mergers, i.e. between kicks. The values of 
Ti/2/trnerge are (0,0.01,0.03,0.1,0.3,1); when this ratio is zero, the cumulative distribution 
is a step-function, i.e. the BH would have returned to the origin. 
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The largest value expected for V is ~ 200 km s" 1 in the absence of spin (Gonzalez et 
al. 2006; Sopuerta, Yunes & Laguna 2006). By comparison, escape velocities from giant 
elliptical galaxies (My < —18) are always greater than ~ 400 km s _1 (Merritt et al. 2004). 
Hence the probability of finding a significant displaced supermassive BH in a bright E galaxy 
is very small, unless a merger occurred very recently. Escape velocities from dwarf elliptical 
(dE) galaxies (—12 < My < —18) are < 150 km s" 1 if only the stellar contribution to the 
potential is considered, and < 300 km s _1 if the dark matter potential is included (Merritt 
et al. 2004). (Escape velocities from luminous E galaxies are dominated by the stars.) 
However in these galaxies the mean time between mergers is believed to be very long, again 
implying a low probability of finding a displaced BH. It is, anyway, worth noting that BH 
return times would be lengthened in galaxies containing DM halos. Actually, in luminous 
galaxies, central escape velocities are determined essentially by the stellar distribution; the 
DM halos have only a small effect (Merritt et al. 2004). On the other side, escape velocities 
in dE galaxies can be significantly affected by DM, as mentioned in the final paragraph of 
previous Section 4. At this regard, useful contribution will come by the quantitative study 
of dynamical friction in self-consistent models of triaxial cuspy galaxies with dark matter 
haloes, as those obtained by Capuzzo-Dolcetta et al. (2007) using the Schwarzschild's (1979) 
orbital superpoasition technique. Preliminary results are those of a significant difference in 
the decay times for energies sufficiently high to appreciate the role of the (spatially extended) 
dark matter distribution. 

As argued above, a BH kicked with velocity V will carry with it material that was 
orbiting with velocity v > V before the kick. The size of the region containing this mass is 



with M 8 = M./1O 8 M and 0200 = o"/200 km s" 1 ; we have set v e « 4<r, appropriate for the 
center of a galaxy. A substantial displacement requires V ~ v e , hence the sphere of entrained 
mass will be of order a parsec in radius. This is sufficient to include the inner accretion disk 
and much of the broad-line region gas, implying that a kicked BH can continue shining 
for some time as a "naked" quasar. This interpretation has in fact been advanced for the 
quasar associated with the HE0450-2958 system, which appears to lack a luminous host 
galaxy (Magain et al. 2005; Haehnelt et al. 2005; Hoffman & Loeb 2006). However a 
number of arguments suggest that the kick hypothesis is unlikely (Merritt et al. 2005); for 
instance, the quasar spectrum exhibits lines associated with the narrow emission line region 
at distances too great to have remain bound to an ejected BH. A recent re-analysis of the 



(Eq. 28) 




(36b) 



(36a) 
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quasar image (Kim et al. 2006) also suggests that the presence of a host galaxy can not 
be ruled out. Nevertheless, detection of a broad emission line spectrum from gas that is 
displaced spatially or kinematically from the center of a galaxy would be strong evidence for 
a kick, particularly if the host galaxy exhibited additional signs of a recent merger. 

5. Conclusions 

We integrated the motion of "kicked" BHs in triaxial models of galaxies, using the 
expressions derived by Pesce, Capuzzo-Dolcetta & Vietri (1992) and Capuzzo-Dolcetta & 
Vicari (2006) for the dynamical friction force in an anisotropic stellar distribution. The 
velocity dispersion components were computed from fully self-consistent triaxial models, 
constructed via orbital superposition. We considered different possible forms for the stellar 
density at the center of the galaxy, since ejection of the BH would significantly affect the 
distribution of stars there. Our main results can be summarized as follows: 

1. Dynamical friction increases the effective escape velocity from a galaxy. This effect 
is modest, roughly a few percent, in galaxy models with shallow central density profiles, but 
can be very significant in galaxies with p ~ r~ 2 central density profiles, since the frictional 
force acting on the BH is so strong near the center (Tabled]). 

2. Since the dynamical friction force in a triaxial galaxy depends on angle as well 
as distance from the center, escape velocities are a function of "launching angle", being 
greatest in the direction of the long axis. Again, this effect is modest in models with low 
central concentration but can be appreciable in galaxies with high central densities (Fig. 

3. The time for a kicked BH to return to the center with zero velocity is longer in a 
triaxial galaxy than in a spherical galaxy with the same radial density profile and length 
scale (a&c) 1 / 3 , and when kicked with the same fraction of the escape velocity. The main 
reason is that trajectories in the triaxial geometry are not linear (unless they are exactly 
along the coordinate axes) and a kicked BH does not return precisely through the center, 
thus reducing the average dynamical friction force (Figs. [3] H]). Infall times are typically 
several times longer than in the spherical geometry (Figs. [MS]). 

4. In spite of the delaying effects of the triaxial geometry, BHs with masses similar to 
those observed in real galaxies (M pa 10~ 3 M ga i) return to the center in less than ~ a galaxy 
crossing time, unless the kick velocity is a large fraction of the escape velocity. Since escape 
velocities in giant elliptical galaxies are large compared with the maximum kick velocities 
so far computed by numerical relativists, the chance of finding a BH substantially displaced 
from the center of such a galaxy is small (Figs. 13, 14). Escape velocities are smaller in 
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dE galaxies but the mean time between mergers is probably long, again implying a small 
probability of finding a displaced BH. 
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